LAMMPS (30 Jun 2020)
# Demonstrate bispectrum computes

# initialize simulation

variable 	nsteps index 0
variable 	nrep equal 1
variable 	a equal 2.0
units		metal

# generate the box and atom positions using a BCC lattice

variable 	nx equal ${nrep}
variable 	nx equal 1
variable 	ny equal ${nrep}
variable 	ny equal 1
variable 	nz equal ${nrep}
variable 	nz equal 1

boundary	p p p

atom_modify	map hash
lattice         bcc $a
lattice         bcc 2
Lattice spacing in x,y,z = 2 2 2
region		box block 0 ${nx} 0 ${ny} 0 ${nz}
region		box block 0 1 0 ${ny} 0 ${nz}
region		box block 0 1 0 1 0 ${nz}
region		box block 0 1 0 1 0 1
create_box	2 box
Created orthogonal box = (0 0 0) to (2 2 2)
  1 by 1 by 1 MPI processor grid
create_atoms	2 box
Created 2 atoms
  create_atoms CPU = 0.000 seconds

mass 		* 180.88

displace_atoms 	all random 0.1 0.1 0.1 123456

# set up reference potential

variable 	zblcutinner equal 4
variable 	zblcutouter equal 4.8
variable 	zblz equal 73
pair_style 	zbl ${zblcutinner} ${zblcutouter}
pair_style 	zbl 4 ${zblcutouter}
pair_style 	zbl 4 4.8
pair_coeff 	* * ${zblz} ${zblz}
pair_coeff 	* * 73 ${zblz}
pair_coeff 	* * 73 73

# choose SNA parameters

variable 	twojmax equal 2
variable 	rcutfac equal 1.0
variable 	rfac0 equal 0.99363
variable 	rmin0 equal 0
variable 	radelem1 equal 2.3
variable 	radelem2 equal 2.0
variable	wj1 equal 1.0
variable	wj2 equal 0.96
variable	quadratic equal 0
variable	bzero equal 0
variable	switch equal 0
variable 	snap_options string "${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}"
1 ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 2.3 ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 2.3 2 ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 2.3 2 1 ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 2.3 2 1 0.96 rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag ${bzero} switchflag ${switch}
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag ${switch}
1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0

# set up per-atom computes

compute 	b all sna/atom ${snap_options}
compute 	b all sna/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
compute 	vb all snav/atom ${snap_options}
compute 	vb all snav/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0
compute 	db all snad/atom ${snap_options}
compute 	db all snad/atom 1 0.99363 2 2.3 2 1 0.96 rmin0 0 quadraticflag 0 bzeroflag 0 switchflag 0

# perform sums over atoms

group 		snapgroup1 type 1
0 atoms in group snapgroup1
group 		snapgroup2 type 2
2 atoms in group snapgroup2
compute         bsum1 snapgroup1 reduce sum c_b[*]
compute         bsum2 snapgroup2 reduce sum c_b[*]
# fix 		bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector
# fix 		bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector
compute		vbsum all reduce sum c_vb[*]
# fix 		vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector
variable	db_2_25 equal c_db[2][25]

thermo 		100

# test output:   1: total potential energy
#                2: xy component of stress tensor
#                3: Sum(B_{000}^i, all i of type 2)
#                4: xz component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j)
#                5: y component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2)
#
#                followed by 5 counterparts from compute snap

# run compute mliap with gradgradflag = 1

compute   	snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 1
SNAP keyword rcutfac 1.0 
SNAP keyword twojmax 2 
SNAP keyword nelems 2 
SNAP keyword elems A 
SNAP keyword radelems 2.3 
SNAP keyword welems 1.0 
SNAP keyword rfac0 0.99363 
SNAP keyword rmin0 0 
SNAP keyword bzeroflag 0 
SNAP keyword switchflag 0 
fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.gg1.dat mode vector

thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 		c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
thermo_modify 	norm no

run             ${nsteps}
run             0
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 6.8
  ghost atom cutoff = 6.8
  binsize = 3.4, bins = 1 1 1
  5 neighbor lists, perpetual/occasional/extra = 1 4 0
  (1) pair zbl, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d/newton
      bin: standard
  (2) compute sna/atom, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
  (3) compute snav/atom, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
  (4) compute snad/atom, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
  (5) compute mliap, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 10.75 | 10.75 | 10.75 Mbytes
PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 
   322.86952    1505558.1    364182.88   -240.25066    1381.7961    322.86952    1505558.1    364182.88   -240.25066    1381.7961 
Loop time of 0 on 1 procs for 0 steps with 2 atoms

0.0% CPU use with 1 MPI tasks x no OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0          | 0          | 0          |   0.0 |  0.00
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 0          | 0          | 0          |   0.0 |  0.00
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 0          |            |       |  0.00

Nlocal:    2.0 ave 2.0 max 2.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    853.0 ave 853.0 max 853.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    330.0 ave 330.0 max 330.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs:  660.0 ave 660.0 max 660.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 660
Ave neighs/atom = 330.0
Neighbor list builds = 0
Dangerous builds = 0

uncompute	snap
unfix		snap

# run compute mliap with gradgradflag = 0

compute   	snap all mliap descriptor sna compute.mliap.descriptor model linear gradgradflag 0
SNAP keyword rcutfac 1.0 
SNAP keyword twojmax 2 
SNAP keyword nelems 2 
SNAP keyword elems A 
SNAP keyword radelems 2.3 
SNAP keyword welems 1.0 
SNAP keyword rfac0 0.99363 
SNAP keyword rmin0 0 
SNAP keyword bzeroflag 0 
SNAP keyword switchflag 0 
fix 		snap all ave/time 1 1 1 c_snap[*] file compute.snap.gg0.dat mode vector

thermo_style	custom 		pe            pxy            c_bsum2[1]   c_vbsum[55]    v_db_2_25 		c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12]
WARNING: New thermo_style command, previous thermo_modify settings will be lost (../output.cpp:694)
thermo_modify 	norm no

run             ${nsteps}
run             0
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 6.8
  ghost atom cutoff = 6.8
  binsize = 3.4, bins = 1 1 1
  5 neighbor lists, perpetual/occasional/extra = 1 4 0
  (1) pair zbl, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d/newton
      bin: standard
  (2) compute sna/atom, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
  (3) compute snav/atom, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
  (4) compute snad/atom, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
  (5) compute mliap, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 22.89 | 22.89 | 22.89 Mbytes
PotEng Pxy c_bsum2[1] c_vbsum[55] v_db_2_25 c_snap[1][13] c_snap[13][13] c_snap[1][8] c_snap[12][12] c_snap[6][12] 
   322.86952    1505558.1    364182.88   -240.25066    1381.7961    322.86952    1505558.1    364182.88   -240.25066    1381.7961 
Loop time of 1e-06 on 1 procs for 0 steps with 2 atoms

100.0% CPU use with 1 MPI tasks x no OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0          | 0          | 0          |   0.0 |  0.00
Neigh   | 0          | 0          | 0          |   0.0 |  0.00
Comm    | 0          | 0          | 0          |   0.0 |  0.00
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 1e-06      |            |       |100.00

Nlocal:    2.0 ave 2.0 max 2.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    853.0 ave 853.0 max 853.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    330.0 ave 330.0 max 330.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs:  660.0 ave 660.0 max 660.0 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 660
Ave neighs/atom = 330.0
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00
